魔数0x5f3759df - [2026重制版]

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );        // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
    // y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, can be removed

    return y;
}

┌─── Sign ───┬──── Exponent ────┬───────── Mantissa ─────────┐
│    1 bit   │     8 bits       │          23 bits            │
│     S      │        E         │             M               │
└────────────┴──────────────────┴─────────────────────────────┘

#include <immintrin.h>
#include <cmath>

// SSE/AVX版本的快速平方根倒数
__m256 fast_rsqrt_ps(__m256 x) {
    // 使用AVX的rsqrt指令（硬件支持，1个周期）
    __m256 approx = _mm256_rsqrt_ps(x);

    // 牛顿迭代精化（可选，提高精度）
    __m256 x2 = _mm256_mul_ps(x, _mm256_set1_ps(0.5f));
    __m256 threehalfs = _mm256_set1_ps(1.5f);
    __m256 y2 = _mm256_mul_ps(approx, approx);
    __m256 term = _mm256_sub_ps(threehalfs, _mm256_mul_ps(x2, y2));
    __m256 result = _mm256_mul_ps(approx, term);

    return result;
}

// 使用示例：批量归一化向量
void normalize_vectors_sse(const float* input,
                           float* output,
                           size_t count) {
    size_t simd_end = count - (count % 8);  // AVX一次处理8个float

    for (size_t i = 0; i < simd_end; i += 8) {
        __m256 v = _mm256_loadu_ps(input + i);

        // 计算点积 v·v 并做水平求和
        __m256 sq = _mm256_mul_ps(v, v);
        __m128 lo = _mm256_castps256_ps128(sq);
        __m128 hi = _mm256_extractf128_ps(sq, 1);
        __m128 sum128 = _mm_add_ps(lo, hi);
        __m128 shuf = _mm_shuffle_ps(sum128, sum128, _MM_SHUFFLE(2, 3, 0, 1));
        __m128 sums = _mm_add_ps(sum128, shuf);
        shuf = _mm_shuffle_ps(sums, sums, _MM_SHUFFLE(1, 0, 3, 2));
        sums = _mm_add_ps(sums, shuf);
        float norm_sq = _mm_cvtss_f32(sums);

        // 广播同一个 |v|^2，保证每个元素乘以同一个 1/|v|
        __m256 inv_sqrt = fast_rsqrt_ps(_mm256_set1_ps(norm_sq));

        // 归一化
        __m256 normalized = _mm256_mul_ps(v, inv_sqrt);

        _mm256_storeu_ps(output + i, normalized);
    }

    // 处理剩余元素（标量版本：按尾部元素的平方和整体归一化）
    float tail_sum = 0.0f;
    for (size_t i = simd_end; i < count; i++) {
        tail_sum += input[i] * input[i];
    }
    float inv_len = Q_rsqrt(tail_sum);  // 使用经典的快速算法
    for (size_t i = simd_end; i < count; i++) {
        output[i] = input[i] * inv_len;
    }
}

// HLSL (DirectX) 着色器中的使用
// 现代GPU都有原生rsq (reciprocal square root)指令

cbuffer Constants : register(b0) {
    float4x4 viewMatrix;
    float4 lightPositions[MAX_LIGHTS];
    float4 lightColors[MAX_LIGHTS];
};

struct PSInput {
    float4 position : SV_POSITION;
    float3 normal : NORMAL;
    float3 worldPos : WORLD_POS;
};

float4 main(PSInput input) : SV_TARGET {
    float3 normal = normalize(input.normal);
    float3 color = float3(0.0, 0.0, 0.0);

    for (int i = 0; i < MAX_LIGHTS; i++) {
        float3 lightDir = normalize(lightPositions[i].xyz - input.worldPos);

        // 使用内置的rsqrt进行快速计算
        // 编译后会映射到GPU的RSQ指令（通常1个时钟周期）
        float distance = length(lightPositions[i].xyz - input.worldPos);
        float attenuation = 1.0 / (distance * distance);  // 这里隐含使用了rsqrt

        float diffuse = max(dot(normal, lightDir), 0.0);
        color += diffuse * lightColors[i].rgb * attenuation;
    }

    return float4(color, 1.0);
}

// 如果需要手动实现（某些旧硬件或教学目的）:
float fastRsqrt(float number) {
    // 在HLSL中可以通过asfloat/asint进行类型转换
    int i = asint(number);
    i = 0x5f3759df - (i >> 1);
    float y = asfloat(i);
    y = y * (1.5 - 0.5 * number * y * y);  // 牛顿迭代
    return y;
}

import math
import torch
import torch.nn.functional as F

class FastInverseSqrt(torch.autograd.Function):
    """
    PyTorch自定义算子：快速平方根倒数
    用于推理时的性能优化
    """
    @staticmethod
    def forward(ctx, x):
        # 使用C++/CUDA扩展实现快速算法
        # 这里展示Python等价逻辑
        output = torch.empty_like(x)

        # 对于GPU张量，使用优化的CUDA kernel
        if x.is_cuda:
            import fast_inverse_sqrt_cuda
            output = fast_inverse_sqrt_cuda.forward(x)
        else:
            # CPU fallback：使用PyTorch内置的rsqrt（已高度优化）
            output = torch.rsqrt(x)

        ctx.save_for_backward(x, output)
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, y = ctx.saved_tensors
        # dy/dx = -0.5 * x^(-3/2) = -0.5 * y^3 / x
        # 但更稳定的数值实现：
        grad_input = grad_output * (-0.5 * y * y * y)
        return grad_input


# 在Transformer注意力机制中的应用
class EfficientAttention(torch.nn.Module):
    def __init__(self, d_model, n_heads):
        super().__init__()
        self.d_model = d_model
        self.n_heads = n_heads
        self.d_k = d_model // n_heads

        self.q_proj = torch.nn.Linear(d_model, d_model)
        self.k_proj = torch.nn.Linear(d_model, d_model)
        self.v_proj = torch.nn.Linear(d_model, d_model)
        self.out_proj = torch.nn.Linear(d_model, d_model)

    def forward(self, x, mask=None):
        batch_size = x.size(0)

        # 投影到Q, K, V
        q = self.q_proj(x).view(batch_size, -1, self.n_heads, self.d_k).transpose(1, 2)
        k = self.k_proj(x).view(batch_size, -1, self.n_heads, self.d_k).transpose(1, 2)
        v = self.v_proj(x).view(batch_size, -1, self.n_heads, self.d_k).transpose(1, 2)

        # 计算注意力分数：除以 sqrt(d_k)，即乘以 1/sqrt(d_k)
        # 用 FastInverseSqrt 预计算并缓存这个缩放系数
        scale = FastInverseSqrt.apply(
            torch.tensor(float(self.d_k), device=x.device)
        )
        scores = torch.matmul(q, k.transpose(-2, -1)) * scale

        if mask is not None:
            scores = scores.masked_fill(mask == 0, -1e9)

        attn_weights = F.softmax(scores, dim=-1)

        # 应用注意力权重
        context = torch.matmul(attn_weights, v)

        # 合并多头
        context = context.transpose(1, 2).contiguous().view(
            batch_size, -1, self.d_model
        )

        return self.out_proj(context)


# TensorFlow/XLA中的使用
import tensorflow as tf

@tf.function(jit_compile=True)
def layer_norm_fast(x, epsilon=1e-6):
    """
    使用XLA编译优化的Layer Normalization
    XLA会自动将rsqrt转换为硬件最优指令
    """
    mean = tf.reduce_mean(x, axis=-1, keepdims=True)
    variance = tf.reduce_mean(tf.math.square(x - mean), axis=-1, keepdims=True)

    # rsqrt在XLA下会被编译为最优实现
    inv_std = tf.math.rsqrt(variance + epsilon)

    return (x - mean) * inv_std

// ARM Cortex-M4/M7上的NEON优化版本
// 适用于STM32、nRF52、ESP32等平台
#include <arm_neon.h>

void normalize_vector_neon(const float* input,
                            float* output,
                            uint32_t vector_count) {
    uint32_t simd_count = vector_count & ~3U;  // NEON一次处理4个float

    for (uint32_t i = 0; i < simd_count; i += 4) {
        float32x4_t v = vld1q_f32(input + i);

        // 计算平方和
        float32x4_t sq = vmulq_f32(v, v);

        // 向量归约求和
        float32x2_t pair = vpadd_f32(
            vget_low_f32(sq),
            vget_high_f32(sq)
        );
        float32x2_t sum = vpadd_f32(pair, pair);

        // 提取标量并广播
        float sum_scalar = vget_lane_f32(sum, 0);
        float32x4_t sum_vec = vmovq_n_f32(sum_scalar);

        // 快速平方根倒数
        // ARM NEON有vrsqrte指令（估计 reciprocal sqrt）
        float32x4_t estimate = vrsqrteq_f32(sum_vec);

        // Newton-Raphson迭代精化（1-2次）
        float32x4_t estimate_sq = vmulq_f32(estimate, estimate);
        float32x4_t half_sum = vmulq_n_f32(sum_vec, 0.5f);
        float32x4_t three_halves = vmovq_n_f32(1.5f);
        float32x4_t term = vsubq_f32(three_halves, vmulq_f32(half_sum, estimate_sq));

        // 归一化
        float32x4_t normalized = vmulq_f32(v, vmulq_f32(estimate, term));

        vst1q_f32(output + i, normalized);
    }

    // 处理剩余向量（标量版本）
    for (uint32_t i = simd_count; i < vector_count; i++) {
        float x = input[i];
        float y = Q_rsqrt(x * x);  // 使用经典的快速算法
        output[i] = x * y;
    }
}

// ESP-IDF (ESP32) 上的定点数版本
// 用于超低功耗传感器数据处理
#include <esp_timer.h>
#include <math.h>

#define FIXED_POINT_SHIFT  16
#define FLOAT_TO_FIXED(x)  ((int32_t)((x) * (1 << FIXED_POINT_SHIFT)))
#define FIXED_TO_FLOAT(x)  ((float)(x) / (1 << FIXED_POINT_SHIFT))

// 定点数版本的快速平方根倒数
// 避免昂贵的浮点运算，适用于ESP32的FPU未启用时
int32_t fixed_q_rsqrt(int32_t number) {
    // 将定点数转换为伪浮点表示
    int32_t x2 = (number >> 1);  // 除以2
    int32_t y = number;

    // 使用移位模拟浮点的位级hack
    // 注意：以下仅为思路示意，并不能直接得到正确结果；
    // 实际的定点实现需要按对数近似重新推导（可参考 McEniry/Lomont 的推导）
    int32_t i = (y >> 15) | 0x5f400000;  // 调整后的魔数
    y = i - (i >> 1);

    // 定点数牛顿迭代
    y = (y * ((3 << FIXED_POINT_SHIFT) - (x2 * y * y >> FIXED_POINT_SHIFT))) >> (FIXED_POINT_SHIFT + 1);

    return y;
}

// IMU传感器数据归一化（ESP32应用示例）
void normalize_imu_data(int16_t* accel, float* norm_accel) {
    // 加速度计原始数据通常是12-16位整数
    int32_t x_sq = (int32_t)accel[0] * accel[0];
    int32_t y_sq = (int32_t)accel[1] * accel[1];
    int32_t z_sq = (int32_t)accel[2] * accel[2];

    int32_t magnitude_sq = x_sq + y_sq + z_sq;

    // 使用快速算法计算1/sqrt(magnitude_sq)
    int32_t inv_magnitude = fixed_q_rsqrt(magnitude_sq);

    // 归一化
    norm_accel[0] = FIXED_TO_FLOAT(accel[0] * inv_magnitude);
    norm_accel[1] = FIXED_TO_FLOAT(accel[1] * inv_magnitude);
    norm_accel[2] = FIXED_TO_FLOAT(accel[2] * inv_magnitude);
}

"""
2026年各平台的平方根倒数性能对比测试
"""
import timeit
import numpy as np
import math

# 方法1：Python标准库math.sqrt
def method_python_math(x):
    return 1.0 / math.sqrt(x)

# 方法2：NumPy向量化
def method_numpy(x):
    return 1.0 / np.sqrt(np.array([x]))[0]

# 方法3：快速近似算法（纯Python模拟）
def method_fast_approx(x):
    # 模拟C语言的位操作
    import struct
    x2 = x * 0.5
    y = x
    packed = struct.pack('!f', y)
    i = struct.unpack('!i', packed)[0]
    i = 0x5f3759df - (i >> 1)
    packed = struct.pack('!i', i)
    y = struct.unpack('!f', packed)[0]
    y = y * (1.5 - x2 * y * y)
    return y

# 测试函数
def benchmark():
    test_values = [0.5, 1.0, 2.0, 5.0, 10.0, 100.0, 1234.5678]

    print("=== 平方根倒数算法性能对比 ===\n")
    print(f"{'方法':<25} {'时间(μs)':<12} {'相对误差':<12}")
    print("-" * 50)

    for name, func in [
        ("math.sqrt", method_python_math),
        ("numpy.sqrt", method_numpy),
        ("Fast Approx", method_fast_approx),
    ]:
        times = []
        errors = []

        for val in test_values:
            # 时间测量
            t = timeit.timeit(lambda: func(val), number=10000)
            times.append(t * 100)  # 转换为微秒

            # 精度测量
            exact = 1.0 / math.sqrt(val)
            approx = func(val)
            rel_error = abs(exact - approx) / exact * 100
            errors.append(rel_error)

        avg_time = np.mean(times)
        avg_error = np.mean(errors)
        print(f"{name:<25} {avg_time:<12.4f} {avg_error:<12.6f}%")

    print("\n=== 各平台实际性能（理论值）===")
    print("""
    平台              | 延迟(ns) | 吞吐(GOPS) | 适用场景
    ------------------|-----------|------------|------------------
    Intel AVX-512     | ~0.5      | ~2000      | HPC, 数据中心
    NVIDIA A100 GPU   | ~0.1      | ~20000     | AI训练/推理
    ARM Cortex-M7 FPU | ~20       | ~50        | 嵌入式系统
    ESP32 (软件模拟)  | ~500      | ~2         | IoT传感器节点
    Quake III时代CPU  | ~100      | ~10        | 历史参考
    """)


if __name__ == "__main__":
    benchmark()